"""
    FlowDependentCorrelations

Sophisticated flow-dependent error correlation modeling for atmospheric data 
assimilation in GSI. This module implements advanced methods for adapting 
background error correlations based on atmospheric flow characteristics,
providing physically consistent error structures that evolve with meteorological
conditions.

Flow-dependent correlations address fundamental limitations of static climatological
error covariances by:
- Adapting to synoptic-scale flow patterns and weather regimes
- Following atmospheric transport and mixing processes
- Respecting conservation laws and dynamical constraints
- Incorporating multi-scale flow interactions
- Supporting both geostrophic and ageostrophic flow components

Mathematical Framework:
Flow-dependent correlation functions are modeled as:
```
C_flow(r₁,r₂,t) = C_static(r₁,r₂) · F(ψ,χ,∇·v,ζ) · T(Δt)
```

Where:
- C_static is the climatological background correlation
- F depends on stream function ψ, velocity potential χ, divergence ∇·v, vorticity ζ
- T(Δt) accounts for temporal evolution of correlations
- Flow-dependent modifications preserve correlation structure positivity

Key Flow Characteristics:
- Stream function analysis for rotational flow components
- Velocity potential decomposition for divergent flow
- Deformation analysis for stretching and shearing effects
- Stability measures for convective and baroclinic instability
- Transport pathway tracking for advective correlation evolution

Applications:
- Hurricane and tropical cyclone error modeling
- Frontal zone correlation enhancement
- Convective-scale flow-dependent structures
- Jet stream and upper-level dynamics
- Multi-scale nested domain consistency
"""
module FlowDependentCorrelations

using LinearAlgebra
using SparseArrays
using StaticArrays

export FlowDependentCorrelation, FlowCharacteristics, StreamVorticityAnalysis
export GeostrophicAgeostrophicDecomposition, FlowStabilityAnalysis
export TransportPathwayTracker, FlowRegimeClassification
export compute_flow_dependent_correlations, decompose_velocity_field
export analyze_flow_stability, track_transport_pathways, classify_flow_regime

"""
    FlowRegime

Enumeration of atmospheric flow regimes for correlation adaptation.
"""
@enum FlowRegime begin
    QUIESCENT_FLOW      # Weak gradient, light winds
    GEOSTROPHIC_FLOW    # Balanced geostrophic flow
    AGEOSTROPHIC_FLOW   # Unbalanced, ageostrophic components
    CYCLONIC_FLOW       # Cyclonic circulation patterns
    ANTICYCLONIC_FLOW   # Anticyclonic circulation patterns
    FRONTAL_ZONE        # Sharp gradient/frontal regions
    JET_STREAM          # High-speed jet stream flow
    CONVECTIVE_FLOW     # Convectively unstable conditions
    BAROCLINIC_ZONE     # Baroclinically unstable regions
    BLOCKING_PATTERN    # Atmospheric blocking conditions
end

"""
    FlowCharacteristics{T<:AbstractFloat}

Comprehensive flow field analysis results for correlation modeling.

# Fields
- `stream_function::Matrix{T}`: Stream function ψ (rotational component)
- `velocity_potential::Matrix{T}`: Velocity potential χ (divergent component)
- `relative_vorticity::Matrix{T}`: Relative vorticity ζ = ∂v/∂x - ∂u/∂y
- `divergence::Matrix{T}`: Horizontal divergence ∇·v = ∂u/∂x + ∂v/∂y
- `deformation_magnitude::Matrix{T}`: Total deformation magnitude
- `stretching_deformation::Matrix{T}`: Stretching deformation component
- `shearing_deformation::Matrix{T}`: Shearing deformation component
- `okubo_weiss_parameter::Matrix{T}`: Okubo-Weiss parameter (strain vs vorticity)
- `flow_regime_classification::Matrix{FlowRegime}`: Local flow regime at each point
"""
struct FlowCharacteristics{T<:AbstractFloat}
    stream_function::Matrix{T}
    velocity_potential::Matrix{T}
    relative_vorticity::Matrix{T}
    divergence::Matrix{T}
    deformation_magnitude::Matrix{T}
    stretching_deformation::Matrix{T}
    shearing_deformation::Matrix{T}
    okubo_weiss_parameter::Matrix{T}
    flow_regime_classification::Matrix{FlowRegime}
end

"""
    StreamVorticityAnalysis{T<:AbstractFloat}

Advanced stream function and vorticity analysis for flow decomposition.

# Fields
- `grid_spacing::Tuple{T,T}`: Grid spacing (dx, dy) in meters
- `coriolis_parameter::Matrix{T}`: Coriolis parameter field f(lat)
- `laplacian_operator::SparseMatrixCSC{T,Int}`: 2D Laplacian operator ∇²
- `gradient_operators::Tuple{SparseMatrixCSC{T,Int}, SparseMatrixCSC{T,Int}}`: ∇x, ∇y operators
- `vorticity_solver_method::Symbol`: Solver method (:direct, :iterative, :multigrid)
- `boundary_conditions::Symbol`: Boundary condition treatment
"""
struct StreamVorticityAnalysis{T<:AbstractFloat}
    grid_spacing::Tuple{T,T}
    coriolis_parameter::Matrix{T}
    laplacian_operator::SparseMatrixCSC{T,Int}
    gradient_operators::Tuple{SparseMatrixCSC{T,Int}, SparseMatrixCSC{T,Int}}
    vorticity_solver_method::Symbol
    boundary_conditions::Symbol
    
    function StreamVorticityAnalysis{T}(;
        grid_size::Tuple{Int,Int},
        grid_spacing::Tuple{T,T} = (T(50000), T(50000)),  # 50 km default
        latitude_field::Union{Matrix{T}, Nothing} = nothing,
        vorticity_solver_method::Symbol = :iterative,
        boundary_conditions::Symbol = :zero_gradient
    ) where T
        nx, ny = grid_size
        
        # Create Coriolis parameter field
        if latitude_field === nothing
            # Default latitude field (simplified)
            lats = range(T(-90), T(90), length=ny) .* π/180
            latitude_field = repeat(reshape(lats, 1, :), nx, 1)
        end
        
        omega = T(7.2921e-5)  # Earth rotation rate
        coriolis_f = 2 * omega .* sin.(latitude_field)
        
        # Create finite difference operators
        lapl_op = create_laplacian_operator(T, nx, ny, grid_spacing)
        grad_x, grad_y = create_gradient_operators(T, nx, ny, grid_spacing)
        
        new{T}(grid_spacing, coriolis_f, lapl_op, (grad_x, grad_y), 
               vorticity_solver_method, boundary_conditions)
    end
end

"""
    GeostrophicAgeostrophicDecomposition{T<:AbstractFloat}

Decomposition of velocity field into geostrophic and ageostrophic components.

# Fields
- `geostrophic_u::Matrix{T}`: Geostrophic zonal velocity u_g
- `geostrophic_v::Matrix{T}`: Geostrophic meridional velocity v_g
- `ageostrophic_u::Matrix{T}`: Ageostrophic zonal velocity u_a = u - u_g
- `ageostrophic_v::Matrix{T}`: Ageostrophic meridional velocity v_a = v - v_g
- `geostrophic_streamfunction::Matrix{T}`: Geostrophic stream function ψ_g
- `total_kinetic_energy::Matrix{T}`: Total kinetic energy ½(u² + v²)
- `geostrophic_ke_fraction::Matrix{T}`: Fraction of KE in geostrophic component
"""
struct GeostrophicAgeostrophicDecomposition{T<:AbstractFloat}
    geostrophic_u::Matrix{T}
    geostrophic_v::Matrix{T}
    ageostrophic_u::Matrix{T}
    ageostrophic_v::Matrix{T}
    geostrophic_streamfunction::Matrix{T}
    total_kinetic_energy::Matrix{T}
    geostrophic_ke_fraction::Matrix{T}
end

"""
    FlowStabilityAnalysis{T<:AbstractFloat}

Analysis of flow stability characteristics for correlation modeling.

# Fields
- `richardson_number::Matrix{T}`: Richardson number Ri = N²/S² 
- `ertel_potential_vorticity::Matrix{T}`: Ertel potential vorticity
- `baroclinic_instability_growth::Matrix{T}`: Baroclinic growth rate estimates
- `convective_instability_index::Matrix{T}`: Convective instability measures
- `static_stability::Matrix{T}`: Static stability parameter N²
- `wind_shear_magnitude::Matrix{T}`: Total wind shear magnitude
"""
struct FlowStabilityAnalysis{T<:AbstractFloat}
    richardson_number::Matrix{T}
    ertel_potential_vorticity::Matrix{T}
    baroclinic_instability_growth::Matrix{T}
    convective_instability_index::Matrix{T}
    static_stability::Matrix{T}
    wind_shear_magnitude::Matrix{T}
end

"""
    TransportPathwayTracker{T<:AbstractFloat}

Lagrangian transport pathway tracking for correlation evolution.

# Fields
- `particle_trajectories::Array{T,3}`: Trajectory coordinates (n_particles × n_time × 2)
- `trajectory_correlations::Matrix{T}`: Correlation evolution along trajectories
- `mixing_coefficients::Matrix{T}`: Turbulent mixing strength estimates
- `transport_barriers::Matrix{Bool}`: Location of transport barriers
- `finite_time_lyapunov_exponent::Matrix{T}`: FTLE field for mixing analysis
"""
struct TransportPathwayTracker{T<:AbstractFloat}
    particle_trajectories::Array{T,3}
    trajectory_correlations::Matrix{T}
    mixing_coefficients::Matrix{T}
    transport_barriers::Matrix{Bool}
    finite_time_lyapunov_exponent::Matrix{T}
end

"""
    FlowDependentCorrelation{T<:AbstractFloat}

Main flow-dependent correlation modeling system.

# Fields
- `stream_vorticity_analyzer::StreamVorticityAnalysis{T}`: Flow field analyzer
- `base_correlations::Matrix{T}`: Static/climatological correlations
- `flow_characteristics::Union{FlowCharacteristics{T}, Nothing}`: Current flow analysis
- `geostrophic_decomposition::Union{GeostrophicAgeostrophicDecomposition{T}, Nothing}`: Velocity decomposition
- `stability_analysis::Union{FlowStabilityAnalysis{T}, Nothing}`: Stability characteristics
- `transport_tracker::Union{TransportPathwayTracker{T}, Nothing}`: Transport pathway tracking
- `flow_adaptation_parameters::Dict{Symbol, T}`: Tuning parameters for flow adaptation
"""
mutable struct FlowDependentCorrelation{T<:AbstractFloat}
    stream_vorticity_analyzer::StreamVorticityAnalysis{T}
    base_correlations::Matrix{T}
    flow_characteristics::Union{FlowCharacteristics{T}, Nothing}
    geostrophic_decomposition::Union{GeostrophicAgeostrophicDecomposition{T}, Nothing}
    stability_analysis::Union{FlowStabilityAnalysis{T}, Nothing}
    transport_tracker::Union{TransportPathwayTracker{T}, Nothing}
    flow_adaptation_parameters::Dict{Symbol, T}
    
    function FlowDependentCorrelation{T}(analyzer::StreamVorticityAnalysis{T}, 
                                        base_correlations::Matrix{T}) where T
        # Default flow adaptation parameters
        flow_params = Dict{Symbol, T}(
            :geostrophic_weight => T(0.8),         # Weight for geostrophic flow adaptation
            :ageostrophic_weight => T(0.6),        # Weight for ageostrophic flow adaptation  
            :vorticity_influence => T(0.5),        # Vorticity influence on correlations
            :divergence_influence => T(0.7),       # Divergence influence on correlations
            :deformation_influence => T(0.4),      # Deformation influence on correlations
            :stability_modulation => T(0.3),       # Stability-based modulation strength
            :transport_decay_rate => T(0.1),       # Transport-based correlation decay
            :regime_adaptation_strength => T(0.6)   # Flow regime adaptation strength
        )
        
        new{T}(analyzer, base_correlations, nothing, nothing, nothing, nothing, flow_params)
    end
end

"""
    decompose_velocity_field(analyzer::StreamVorticityAnalysis{T}, 
                            u_velocity::Matrix{T}, v_velocity::Matrix{T}) where T

Decompose velocity field into rotational and divergent components.

# Arguments
- `analyzer::StreamVorticityAnalysis{T}`: Flow analysis system
- `u_velocity::Matrix{T}`: Zonal velocity component
- `v_velocity::Matrix{T}`: Meridional velocity component

# Returns
- `FlowCharacteristics{T}`: Complete flow field analysis

# Algorithm
1. Compute relative vorticity: ζ = ∂v/∂x - ∂u/∂y
2. Compute divergence: δ = ∂u/∂x + ∂v/∂y  
3. Solve for stream function: ∇²ψ = ζ
4. Solve for velocity potential: ∇²χ = δ
5. Compute deformation tensor and derived quantities
6. Classify flow regime at each grid point
"""
function decompose_velocity_field(analyzer::StreamVorticityAnalysis{T}, 
                                 u_velocity::Matrix{T}, v_velocity::Matrix{T}) where T
    nx, ny = size(u_velocity)
    grad_x, grad_y = analyzer.gradient_operators
    
    # Flatten velocity fields for gradient computation
    u_flat = vec(u_velocity)
    v_flat = vec(v_velocity)
    
    # Compute velocity gradients
    dudx = reshape(grad_x * u_flat, nx, ny)
    dudy = reshape(grad_y * u_flat, nx, ny)
    dvdx = reshape(grad_x * v_flat, nx, ny)
    dvdy = reshape(grad_y * v_flat, nx, ny)
    
    # Compute vorticity and divergence
    relative_vorticity = dvdx - dudy
    divergence = dudx + dvdy
    
    # Solve for stream function: ∇²ψ = ζ
    vorticity_flat = vec(relative_vorticity)
    if analyzer.vorticity_solver_method == :direct
        psi_flat = analyzer.laplacian_operator \ vorticity_flat
    else
        # Iterative solver for large systems
        psi_flat = solve_iterative(analyzer.laplacian_operator, vorticity_flat)
    end
    stream_function = reshape(psi_flat, nx, ny)
    
    # Solve for velocity potential: ∇²χ = δ
    divergence_flat = vec(divergence)
    if analyzer.vorticity_solver_method == :direct
        chi_flat = analyzer.laplacian_operator \ divergence_flat
    else
        chi_flat = solve_iterative(analyzer.laplacian_operator, divergence_flat)
    end
    velocity_potential = reshape(chi_flat, nx, ny)
    
    # Compute deformation tensor components
    # Stretching deformation: def1 = ∂u/∂x - ∂v/∂y
    stretching_deformation = dudx - dvdy
    
    # Shearing deformation: def2 = ∂v/∂x + ∂u/∂y  
    shearing_deformation = dvdx + dudy
    
    # Total deformation magnitude
    deformation_magnitude = sqrt.(stretching_deformation.^2 + shearing_deformation.^2)
    
    # Compute Okubo-Weiss parameter: OW = S² - ζ²
    # where S² = (∂u/∂x - ∂v/∂y)² + (∂v/∂x + ∂u/∂y)²
    strain_squared = stretching_deformation.^2 + shearing_deformation.^2
    okubo_weiss = strain_squared - relative_vorticity.^2
    
    # Classify flow regime at each grid point
    flow_regime_classification = classify_flow_regime_field(u_velocity, v_velocity, 
                                                          relative_vorticity, divergence,
                                                          deformation_magnitude, okubo_weiss)
    
    return FlowCharacteristics{T}(
        stream_function, velocity_potential, relative_vorticity, divergence,
        deformation_magnitude, stretching_deformation, shearing_deformation,
        okubo_weiss, flow_regime_classification
    )
end

"""
    compute_geostrophic_ageostrophic_decomposition(analyzer::StreamVorticityAnalysis{T},
                                                  u_velocity::Matrix{T}, v_velocity::Matrix{T},
                                                  geopotential::Matrix{T}) where T

Decompose velocity into geostrophic and ageostrophic components.
"""
function compute_geostrophic_ageostrophic_decomposition(analyzer::StreamVorticityAnalysis{T},
                                                       u_velocity::Matrix{T}, v_velocity::Matrix{T},
                                                       geopotential::Matrix{T}) where T
    nx, ny = size(u_velocity)
    grad_x, grad_y = analyzer.gradient_operators
    
    # Compute geostrophic winds from geopotential
    # u_g = -(1/f) ∂Φ/∂y,  v_g = (1/f) ∂Φ/∂x
    phi_flat = vec(geopotential)
    dphidx = reshape(grad_x * phi_flat, nx, ny)
    dphidy = reshape(grad_y * phi_flat, nx, ny)
    
    # Geostrophic velocity components
    geostrophic_u = -dphidy ./ analyzer.coriolis_parameter
    geostrophic_v = dphidx ./ analyzer.coriolis_parameter
    
    # Handle near-equatorial regions where f ≈ 0
    small_f_mask = abs.(analyzer.coriolis_parameter) .< 1e-10
    geostrophic_u[small_f_mask] .= u_velocity[small_f_mask]
    geostrophic_v[small_f_mask] .= v_velocity[small_f_mask]
    
    # Ageostrophic components
    ageostrophic_u = u_velocity - geostrophic_u
    ageostrophic_v = v_velocity - geostrophic_v
    
    # Geostrophic stream function
    geostrophic_streamfunction = geopotential ./ analyzer.coriolis_parameter
    geostrophic_streamfunction[small_f_mask] .= T(0)
    
    # Kinetic energy analysis
    total_ke = 0.5 .* (u_velocity.^2 + v_velocity.^2)
    geostrophic_ke = 0.5 .* (geostrophic_u.^2 + geostrophic_v.^2)
    geostrophic_ke_fraction = geostrophic_ke ./ (total_ke .+ T(1e-12))  # Avoid division by zero
    
    return GeostrophicAgeostrophicDecomposition{T}(
        geostrophic_u, geostrophic_v, ageostrophic_u, ageostrophic_v,
        geostrophic_streamfunction, total_ke, geostrophic_ke_fraction
    )
end

"""
    analyze_flow_stability(analyzer::StreamVorticityAnalysis{T}, 
                          u_velocity::Matrix{T}, v_velocity::Matrix{T},
                          temperature::Matrix{T}, pressure_levels::Vector{T}) where T

Analyze flow stability characteristics.
"""
function analyze_flow_stability(analyzer::StreamVorticityAnalysis{T}, 
                               u_velocity::Matrix{T}, v_velocity::Matrix{T},
                               temperature::Matrix{T}, pressure_levels::Vector{T}) where T
    nx, ny = size(u_velocity)
    nz = length(pressure_levels)
    
    # Initialize stability fields
    richardson_number = ones(T, nx, ny) .* T(10)  # Default high Ri (stable)
    ertel_pv = zeros(T, nx, ny)
    baroclinic_growth = zeros(T, nx, ny)
    convective_index = zeros(T, nx, ny)
    static_stability = ones(T, nx, ny) .* T(1e-4)  # Default weak stability
    wind_shear = sqrt.((analyzer.gradient_operators[1] * vec(u_velocity)).^2 + 
                      (analyzer.gradient_operators[2] * vec(v_velocity)).^2)
    wind_shear_magnitude = reshape(wind_shear, nx, ny)
    
    # Simplified stability analysis (would be more complex in full implementation)
    # Richardson number approximation using horizontal shear
    for j in 1:ny, i in 1:nx
        if i > 1 && i < nx && j > 1 && j < ny
            # Horizontal wind shear
            shear_u = (u_velocity[i+1,j] - u_velocity[i-1,j]) / (2 * analyzer.grid_spacing[1])
            shear_v = (v_velocity[i,j+1] - v_velocity[i,j-1]) / (2 * analyzer.grid_spacing[2])
            total_shear = sqrt(shear_u^2 + shear_v^2)
            
            if total_shear > T(1e-6)
                # Simplified Richardson number (would need 3D thermal structure)
                richardson_number[i,j] = static_stability[i,j] / total_shear^2
            end
        end
        
        # Simplified convective index based on local temperature gradients
        if i > 1 && i < nx && j > 1 && j < ny
            temp_variance = var([temperature[i-1,j], temperature[i+1,j], 
                               temperature[i,j-1], temperature[i,j+1]])
            convective_index[i,j] = temp_variance / (temperature[i,j] + T(273.15))
        end
    end
    
    return FlowStabilityAnalysis{T}(
        richardson_number, ertel_pv, baroclinic_growth, 
        convective_index, static_stability, wind_shear_magnitude
    )
end

"""
    compute_flow_dependent_correlations(flow_corr::FlowDependentCorrelation{T},
                                       u_velocity::Matrix{T}, v_velocity::Matrix{T}) where T

Compute flow-dependent correlation matrix modifications.
"""
function compute_flow_dependent_correlations(flow_corr::FlowDependentCorrelation{T},
                                            u_velocity::Matrix{T}, v_velocity::Matrix{T}) where T
    # Analyze flow characteristics
    flow_characteristics = decompose_velocity_field(flow_corr.stream_vorticity_analyzer, 
                                                   u_velocity, v_velocity)
    flow_corr.flow_characteristics = flow_characteristics
    
    nx, ny = size(u_velocity)
    n_points = nx * ny
    
    # Create flow-dependent correlation modifications
    flow_modification = ones(T, n_points, n_points)
    
    # Extract flow parameters
    vort_influence = flow_corr.flow_adaptation_parameters[:vorticity_influence]
    div_influence = flow_corr.flow_adaptation_parameters[:divergence_influence]
    def_influence = flow_corr.flow_adaptation_parameters[:deformation_influence]
    regime_strength = flow_corr.flow_adaptation_parameters[:regime_adaptation_strength]
    
    # Apply flow-dependent modifications
    for j in 1:ny, i in 1:nx
        point_idx = (j-1)*nx + i
        
        # Local flow characteristics at point (i,j)
        local_vorticity = flow_characteristics.relative_vorticity[i,j]
        local_divergence = flow_characteristics.divergence[i,j]
        local_deformation = flow_characteristics.deformation_magnitude[i,j]
        local_regime = flow_characteristics.flow_regime_classification[i,j]
        
        # Compute flow-dependent correlation modifications for surrounding points
        for jj in 1:ny, ii in 1:nx
            target_idx = (jj-1)*nx + ii
            
            if point_idx != target_idx
                # Distance between points
                dx = T(ii - i) * flow_corr.stream_vorticity_analyzer.grid_spacing[1]
                dy = T(jj - j) * flow_corr.stream_vorticity_analyzer.grid_spacing[2]
                distance = sqrt(dx^2 + dy^2)
                
                # Flow-dependent correlation modification
                vorticity_factor = 1 + vort_influence * tanh(local_vorticity * distance / T(100000))
                divergence_factor = 1 + div_influence * tanh(local_divergence * distance / T(100000))
                deformation_factor = 1 + def_influence * tanh(local_deformation * distance / T(100000))
                
                # Regime-specific modifications
                regime_factor = compute_regime_correlation_factor(local_regime, distance, regime_strength)
                
                # Combined modification
                total_modification = vorticity_factor * divergence_factor * 
                                   deformation_factor * regime_factor
                
                # Ensure modification stays within reasonable bounds
                flow_modification[point_idx, target_idx] = clamp(total_modification, T(0.1), T(2.0))
            end
        end
    end
    
    # Apply modifications to base correlations
    modified_correlations = flow_corr.base_correlations .* flow_modification
    
    return modified_correlations
end

"""
    classify_flow_regime_field(u_velocity::Matrix{T}, v_velocity::Matrix{T},
                              vorticity::Matrix{T}, divergence::Matrix{T},
                              deformation::Matrix{T}, okubo_weiss::Matrix{T}) where T

Classify flow regime at each grid point based on flow characteristics.
"""
function classify_flow_regime_field(u_velocity::Matrix{T}, v_velocity::Matrix{T},
                                   vorticity::Matrix{T}, divergence::Matrix{T},
                                   deformation::Matrix{T}, okubo_weiss::Matrix{T}) where T
    nx, ny = size(u_velocity)
    classification = Matrix{FlowRegime}(undef, nx, ny)
    
    # Flow speed for classification
    wind_speed = sqrt.(u_velocity.^2 + v_velocity.^2)
    
    # Classification thresholds (would be tuned based on experience)
    weak_wind_threshold = T(5.0)      # m/s
    strong_wind_threshold = T(25.0)   # m/s
    high_vorticity_threshold = T(1e-4) # 1/s
    high_divergence_threshold = T(1e-4) # 1/s
    high_deformation_threshold = T(1e-4) # 1/s
    
    for j in 1:ny, i in 1:nx
        speed = wind_speed[i,j]
        vort = abs(vorticity[i,j])
        div = abs(divergence[i,j])
        def = deformation[i,j]
        ow = okubo_weiss[i,j]
        
        if speed < weak_wind_threshold
            classification[i,j] = QUIESCENT_FLOW
        elseif speed > strong_wind_threshold
            classification[i,j] = JET_STREAM
        elseif vort > high_vorticity_threshold
            if vorticity[i,j] > 0  # Positive vorticity (cyclonic in NH)
                classification[i,j] = CYCLONIC_FLOW
            else
                classification[i,j] = ANTICYCLONIC_FLOW
            end
        elseif div > high_divergence_threshold
            classification[i,j] = CONVECTIVE_FLOW
        elseif def > high_deformation_threshold
            classification[i,j] = FRONTAL_ZONE
        elseif ow > 0  # Strain-dominated
            classification[i,j] = FRONTAL_ZONE
        elseif ow < -high_vorticity_threshold^2  # Vorticity-dominated
            classification[i,j] = CYCLONIC_FLOW
        else
            classification[i,j] = GEOSTROPHIC_FLOW
        end
    end
    
    return classification
end

"""
    compute_regime_correlation_factor(regime::FlowRegime, distance::T, 
                                    adaptation_strength::T) where T

Compute regime-specific correlation modification factor.
"""
function compute_regime_correlation_factor(regime::FlowRegime, distance::T, 
                                         adaptation_strength::T) where T
    base_factor = T(1)
    
    if regime == QUIESCENT_FLOW
        # Increase correlations slightly for weak gradient flows
        base_factor = 1 + adaptation_strength * T(0.1)
    elseif regime == JET_STREAM
        # Increase along-stream correlations, decrease cross-stream
        base_factor = 1 + adaptation_strength * T(0.2)
    elseif regime == CYCLONIC_FLOW
        # Enhance correlations within vortical structures
        decay_scale = T(200000)  # 200 km
        base_factor = 1 + adaptation_strength * T(0.3) * exp(-distance/decay_scale)
    elseif regime == FRONTAL_ZONE
        # Sharp correlation transitions across fronts
        front_scale = T(50000)   # 50 km
        base_factor = 1 - adaptation_strength * T(0.4) * exp(-distance/front_scale)
    elseif regime == CONVECTIVE_FLOW
        # Localized correlations for convective regimes
        conv_scale = T(20000)    # 20 km
        base_factor = 1 + adaptation_strength * T(0.5) * exp(-distance/conv_scale)
    else
        # Default geostrophic behavior
        base_factor = T(1)
    end
    
    return base_factor
end

"""
    solve_iterative(A::SparseMatrixCSC{T,Int}, b::Vector{T}; 
                   tol::T = T(1e-6), maxiter::Int = 1000) where T

Iterative solver for large sparse linear systems.
"""
function solve_iterative(A::SparseMatrixCSC{T,Int}, b::Vector{T}; 
                        tol::T = T(1e-6), maxiter::Int = 1000) where T
    n = length(b)
    x = zeros(T, n)
    r = b - A * x
    
    # Simple conjugate gradient iteration (would use more sophisticated methods)
    for iter in 1:maxiter
        if norm(r) < tol
            break
        end
        
        # Simplified iteration (would implement full CG)
        alpha = dot(r, r) / dot(r, A * r + T(1e-12) * r)  # Regularization
        x += alpha * r
        r = b - A * x
    end
    
    return x
end

"""
    create_laplacian_operator(::Type{T}, nx::Int, ny::Int, 
                             grid_spacing::Tuple{T,T}) where T

Create 2D Laplacian operator for stream function computation.
"""
function create_laplacian_operator(::Type{T}, nx::Int, ny::Int, 
                                  grid_spacing::Tuple{T,T}) where T
    dx, dy = grid_spacing
    n = nx * ny
    
    I_indices = Int[]
    J_indices = Int[]
    values = T[]
    
    for j in 1:ny, i in 1:nx
        idx = (j-1)*nx + i
        
        # Interior points: 5-point stencil
        if 1 < i < nx && 1 < j < ny
            # Central point coefficient
            push!(I_indices, idx); push!(J_indices, idx); 
            push!(values, -2/dx^2 - 2/dy^2)
            
            # Neighboring points
            push!(I_indices, idx); push!(J_indices, idx-1); push!(values, 1/dx^2)  # Left
            push!(I_indices, idx); push!(J_indices, idx+1); push!(values, 1/dx^2)  # Right
            push!(I_indices, idx); push!(J_indices, idx-nx); push!(values, 1/dy^2) # Bottom
            push!(I_indices, idx); push!(J_indices, idx+nx); push!(values, 1/dy^2) # Top
        else
            # Boundary treatment: identity (would implement proper BCs)
            push!(I_indices, idx); push!(J_indices, idx); push!(values, T(1))
        end
    end
    
    return sparse(I_indices, J_indices, values, n, n)
end

"""
    create_gradient_operators(::Type{T}, nx::Int, ny::Int, 
                             grid_spacing::Tuple{T,T}) where T

Create finite difference gradient operators.
"""
function create_gradient_operators(::Type{T}, nx::Int, ny::Int, 
                                  grid_spacing::Tuple{T,T}) where T
    dx, dy = grid_spacing
    n = nx * ny
    
    # X-gradient (∂/∂x)
    I_x, J_x, V_x = Int[], Int[], T[]
    # Y-gradient (∂/∂y)
    I_y, J_y, V_y = Int[], Int[], T[]
    
    for j in 1:ny, i in 1:nx
        idx = (j-1)*nx + i
        
        # X-gradient
        if 1 < i < nx
            # Central difference
            push!(I_x, idx); push!(J_x, idx-1); push!(V_x, -T(0.5)/dx)
            push!(I_x, idx); push!(J_x, idx+1); push!(V_x, T(0.5)/dx)
        elseif i == 1
            # Forward difference
            push!(I_x, idx); push!(J_x, idx); push!(V_x, -T(1)/dx)
            push!(I_x, idx); push!(J_x, idx+1); push!(V_x, T(1)/dx)
        else # i == nx
            # Backward difference
            push!(I_x, idx); push!(J_x, idx-1); push!(V_x, -T(1)/dx)
            push!(I_x, idx); push!(J_x, idx); push!(V_x, T(1)/dx)
        end
        
        # Y-gradient
        if 1 < j < ny
            # Central difference
            push!(I_y, idx); push!(J_y, idx-nx); push!(V_y, -T(0.5)/dy)
            push!(I_y, idx); push!(J_y, idx+nx); push!(V_y, T(0.5)/dy)
        elseif j == 1
            # Forward difference
            push!(I_y, idx); push!(J_y, idx); push!(V_y, -T(1)/dy)
            push!(I_y, idx); push!(J_y, idx+nx); push!(V_y, T(1)/dy)
        else # j == ny
            # Backward difference
            push!(I_y, idx); push!(J_y, idx-nx); push!(V_y, -T(1)/dy)
            push!(I_y, idx); push!(J_y, idx); push!(V_y, T(1)/dy)
        end
    end
    
    grad_x = sparse(I_x, J_x, V_x, n, n)
    grad_y = sparse(I_y, J_y, V_y, n, n)
    
    return grad_x, grad_y
end

"""
    track_transport_pathways(flow_corr::FlowDependentCorrelation{T},
                           u_velocity::Matrix{T}, v_velocity::Matrix{T}, 
                           time_steps::Vector{T}) where T

Track Lagrangian transport pathways for correlation evolution analysis.
"""
function track_transport_pathways(flow_corr::FlowDependentCorrelation{T},
                                 u_velocity::Matrix{T}, v_velocity::Matrix{T}, 
                                 time_steps::Vector{T}) where T
    nx, ny = size(u_velocity)
    n_time = length(time_steps)
    
    # Initialize particle locations (simplified - would use more sophisticated seeding)
    n_particles = min(100, nx*ny÷10)  # Limit particles for efficiency
    particle_i = rand(1:nx, n_particles)
    particle_j = rand(1:ny, n_particles)
    
    # Track trajectories
    trajectories = zeros(T, n_particles, n_time, 2)
    trajectory_correlations = zeros(T, n_particles, n_time)
    
    # Initialize
    trajectories[:, 1, 1] = T.(particle_i)
    trajectories[:, 1, 2] = T.(particle_j)
    trajectory_correlations[:, 1] .= T(1)  # Perfect correlation at t=0
    
    # Simple Euler integration (would use more sophisticated methods)
    dt = T(3600)  # 1 hour time step
    dx, dy = flow_corr.stream_vorticity_analyzer.grid_spacing
    
    for t in 2:n_time
        for p in 1:n_particles
            # Current particle position
            pi, pj = trajectories[p, t-1, 1], trajectories[p, t-1, 2]
            
            # Bounds checking
            i_idx = max(1, min(nx, round(Int, pi)))
            j_idx = max(1, min(ny, round(Int, pj)))
            
            # Interpolate velocity at particle location
            u_interp = u_velocity[i_idx, j_idx]
            v_interp = v_velocity[i_idx, j_idx]
            
            # Advect particle
            new_pi = pi + u_interp * dt / dx
            new_pj = pj + v_interp * dt / dy
            
            # Store trajectory
            trajectories[p, t, 1] = new_pi
            trajectories[p, t, 2] = new_pj
            
            # Compute correlation decay along trajectory
            # Simplified exponential decay
            decay_rate = flow_corr.flow_adaptation_parameters[:transport_decay_rate]
            trajectory_correlations[p, t] = trajectory_correlations[p, t-1] * 
                                           exp(-decay_rate * dt)
        end
    end
    
    # Create simplified transport tracker
    mixing_coefficients = ones(T, nx, ny) .* T(0.1)  # Simplified
    transport_barriers = falses(nx, ny)              # No barriers for now
    ftle_field = zeros(T, nx, ny)                    # Simplified FTLE
    
    tracker = TransportPathwayTracker{T}(
        trajectories, trajectory_correlations, mixing_coefficients,
        transport_barriers, ftle_field
    )
    
    flow_corr.transport_tracker = tracker
    return tracker
end

"""
    classify_flow_regime(u_velocity::Matrix{T}, v_velocity::Matrix{T}) where T

High-level flow regime classification for entire domain.
"""
function classify_flow_regime(u_velocity::Matrix{T}, v_velocity::Matrix{T}) where T
    # Compute domain-average flow characteristics
    wind_speed = sqrt.(u_velocity.^2 + v_velocity.^2)
    mean_speed = mean(wind_speed)
    max_speed = maximum(wind_speed)
    speed_variability = std(wind_speed)
    
    # Simple domain classification
    if mean_speed < T(3.0)
        return QUIESCENT_FLOW
    elseif max_speed > T(30.0)
        return JET_STREAM
    elseif speed_variability > T(10.0)
        return FRONTAL_ZONE
    else
        return GEOSTROPHIC_FLOW
    end
end

end # module FlowDependentCorrelations